Bayesian Checks

Code

library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)
library(parallel)

theme_set(theme_pubr(legend = "bottom"))

read imputed data

Code

nnns_imputed<- readRDS("../Haojia_work/nnns_imputed.rds")

dat <- lapply(1:20, function(i) complete(nnns_imputed, i))

Use ZOIB package

ZOIB Model

\[ f(y_i|\eta_i) = \begin{cases} p_i, & y_i =0 \\ (1-p_i)q_i, & y_i =1 \\ (1 − p_i)(1 − q_i)Beta(\alpha_{i_1}, \alpha_{i_2}) &y_i \in (0, 1) \end{cases} \]

\[ \begin{aligned} logit\left(\mu^{(0,1)} = \frac{\alpha_1}{\alpha_1+ \alpha_2}\right) = \mathbf x_1 \boldsymbol \beta_1 + I_1(\mathbf z_{1,i}\gamma_{1,i}) \\ log\left(v_i = \alpha_1+ \alpha_2\right) = \mathbf x_2 \boldsymbol \beta_2 + I_2(\mathbf z_{2,i}\gamma_{2,i}) \\ logit\left(p_i\right) = \mathbf x_3 \boldsymbol \beta_3 + I_1(\mathbf z_{3,i}\gamma_{3,i}) \\ logit\left(q_i\right) = \mathbf x_4 \boldsymbol \beta_4 + I_1(\mathbf z_{4,i}\gamma_{4 ,i}) \end{aligned} \]

Fit ZOIB model

Code

set.seed(11282023)

tictoc::tic()

# Define a function for model fitting
fit_model <- function(d) {
  zoib(
    Percent_of_feeds_taken_by_mouth_at_discharge ~
      Pre_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female #x1 design matrix
    | 1 | #x2 design matrix
      Pre_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female | #X3 design matrix
      Pre_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female, #x4 design matrix
    data = d,
    n.response = 1,
    zero.inflation = TRUE,
    one.inflation = TRUE,
    link.mu = "logit",
    link.x0 = "logit",
    link.x1 = "logit",
    random = 0,
    n.chain = 4,
    n.iter = 3000,
    n.thin = 2,
    n.burn = 200,
    seeds = c(11, 29, 20, 23)
  )
}
model_results <- lapply(dat, fit_model)
# Save results
saveRDS(model_results, "pre_op_models.rds")

tictoc::toc()

Code

set.seed(11282023)

tictoc::tic()

# Define a function for model fitting
fit_model <- function(d) {
  zoib(
    Percent_of_feeds_taken_by_mouth_at_discharge ~
      Post_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female #x1 design matrix
    | 1 | #x2 design matrix
      Post_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female | #X3 design matrix
      Post_Op_NNNS_attention_score +
      Length_of_intubation_days +
      Cardiac_Anatomy +
      Age_at_Surgery_days +
      Female, #x4 design matrix
    data = d,
    n.response = 1,
    zero.inflation = TRUE,
    one.inflation = TRUE,
    link.mu = "logit",
    link.x0 = "logit",
    link.x1 = "logit",
    random = 0,
    n.chain = 4,
    n.iter = 3000,
    n.thin = 2,
    n.burn = 200,
    seeds = c(11, 29, 20, 23)
  )
}

# Parallelize model fitting
model_results <- lapply(dat, fit_model)

# Save results
saveRDS(model_results, "post_op_models.rds")

tictoc::toc()

Interpreting output:

b: vector of estimates from Eqn 1; that is, g(mu) = xb*b +z*gamma

d: vector of estimates from Eqn 2; that is, log(eta) = xd*d+z*gamma

b0: vector of estimates from Eqn 3; that is, g(p0) = x0*b0 +z*gamma

b1: vector of estimates from Eqn ;4 that is, g(p1) = x1*b1+z*gamma

Pool chains

Code

post_op_models =readRDS(file="post_op_models.rds")
pre_op_models =readRDS(file="pre_op_models.rds")

post_op_coeff = list()
pre_op_coeff = list()
for(i in 1:length(post_op_models)){
  pre_op_coeff[[i]] = pre_op_models[[i]]$coeff
  post_op_coeff[[i]] = post_op_models[[i]]$coeff
}
pooled_pre_op = runjags::combine.mcmc(mcmc.objects = pre_op_coeff, collapse.chains = FALSE)
pooled_post_op = runjags::combine.mcmc(mcmc.objects = post_op_coeff, collapse.chains = FALSE)

saveRDS(pooled_pre_op, "pooled_pre_op.rds")
saveRDS(pooled_post_op, "pooled_post_op.rds")

Code

pooled_pre_op = readRDS("pooled_pre_op.rds")
pooled_post_op = readRDS("pooled_post_op.rds")

Check convergence

Traceplots

Pre-op

Code

pooled_pre_op |> traceplot()

Post-op

Code

pooled_post_op |> traceplot()

Autocorrelation plots

Pre-op

Code

pooled_pre_op |> autocorr.plot()

Post-op

Code

pooled_post_op |> autocorr.plot()

Summary table

Code

summary_pooled_pre_op = pooled_pre_op |> summary()
summary_pooled_post_op = pooled_post_op |> summary()

rnames = c("Baseline","Attention Score","Length of Intubation (d)",
           "Single Ventricle w/ Arch Obstruction","Two Ventricles w/ Arch Obstruction",
           "Age", "Female")

pre_mean = summary_pooled_pre_op$statistics[,"Mean"]
pre_lb = summary_pooled_pre_op$quantiles[,"2.5%"] 
pre_ub = summary_pooled_pre_op$quantiles[,"97.5%"]



post_mean = summary_pooled_post_op$statistics[,"Mean"]
post_lb   = summary_pooled_post_op$quantiles[,"2.5%"] 
post_ub   = summary_pooled_post_op$quantiles[,"97.5%"] 


pre_b_df= cbind("Mean"= pre_mean[1:7] |> exp() |> round(2),
                 "2.5%"= pre_lb[1:7] |> exp() |> round(2),
                 "97.5%" = pre_ub[1:7] |> exp() |> round(2))

pre_b0_df= cbind("Mean"= pre_mean[8:14] |> exp() |> round(2),
               "2.5%"= pre_lb[8:14] |> exp() |> round(2),
               "97.5%" = pre_ub[8:14] |> exp() |> round(2))
pre_b1_df= cbind("Mean"= pre_mean[15:21] |> exp() |> round(2),
               "2.5%"= pre_lb[15:21] |> exp() |> round(2),
               "97.5%" = pre_ub[15:21] |> exp() |> round(2))

pre_df = cbind(pre_b_df,pre_b0_df,pre_b1_df)
rownames(pre_df) = rnames



post_b_df= cbind("Mean"= post_mean[1:7] |> exp() |> round(2),
                 "2.5%"= post_lb[1:7] |> exp() |> round(2),
                 "97.5%" = post_ub[1:7] |> exp() |> round(2))

post_b0_df= cbind("Mean"= post_mean[8:14] |> exp() |> round(2),
               "2.5%"= post_lb[8:14] |> exp() |> round(2),
               "97.5%" = post_ub[8:14] |> exp() |> round(2))
post_b1_df= cbind("Mean "= post_mean[15:21] |> exp() |> round(2),
               "2.5%"= post_lb[15:21] |> exp() |> round(2),
               "97.5%" = post_ub[15:21] |> exp() |> round(2))

post_df = cbind(post_b_df,post_b0_df,post_b1_df)
rownames(post_df) = rnames


pre_kable = pre_df |>
  kable(caption = "Odds Ratio of Percent Oral Feed Model Results (Pre-Operation)") |>
  add_header_above(header = c("Predictor" = 1,
                            "Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
                              "Odds of 0% Oral Feed" = 3,
                              "Odds of 100% Oral Feed" =3)) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")) |>
  add_footnote(paste("Posterior variance is estimated to be ", pre_mean[22] |> exp() |> round(2),
                     " (",pre_lb[22] |> exp() |> round(2),",", pre_ub[22]|> exp()|> round(2),")"))

post_kable = post_df |>
  kable(caption = "Odds Ratio of Percent Oral Feed Model Results (Post-Operation)") |>
  add_header_above(header = c("Predictor" = 1,
                            "Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
                              "Odds of 0% Oral Feed" = 3,
                              "Odds of 100% Oral Feed" =3)) |>
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))|>
  add_footnote(paste("Posterior variance is estimated to be ", post_mean[22] |> exp() |> round(2),
                     " (",post_lb[22] |> exp() |> round(2),",", post_ub[22]|> exp()|> round(2),")"))

Pre-operation

Code

pre_kable

Percent Oral Feed Model Results (Pre-Operation)

Predictor	Odds of Oral Feed when oral feed between 0 and 1			Odds of 0% Oral Feed			Odds of 100% Oral Feed
	Mean	2.5%	97.5	Mean	2.5%	97.5	Mean	2.5%	97.5
Baseline	0.96	0.21	4.47	0.20	0.02	2.32	0.02	0.00	1.11
Attention Score	1.02	0.71	1.47	0.74	0.37	1.35	1.74	0.63	5.34
Length of Intubation (d)	0.86	0.76	0.98	1.35	1.13	1.67	0.69	0.45	0.99
Single Ventricle w/ Arch Obstruction	1.57	0.71	3.36	6.19	1.90	22.14	0.76	0.03	10.62
Two Ventricles w/ Arch Obstruction	0.81	0.42	1.59	3.43	1.11	11.58	0.52	0.08	3.05
Age	0.98	0.92	1.04	0.97	0.87	1.07	1.19	1.02	1.43
Female	0.65	0.37	1.13	0.56	0.21	1.44	1.94	0.40	9.96

Post-operation

Code

post_kable

Percent Oral Feed Model Results (Post-Operation)

Predictor	Odds of Oral Feed when oral feed between 0 and 1			Odds of 0% Oral Feed			Odds of 100% Oral Feed
	Mean	2.5%	97.5%	Mean	2.5%	97.5%	Mean	2.5%	97.5%
Baseline	0.39	0.04	3.29	0.03	0.00	0.67	0.02	0.00	37.09
Attention Score	1.20	0.82	1.75	1.19	0.71	2.01	1.33	0.34	4.09
Length of Intubation (d)	0.88	0.77	1.00	1.34	1.12	1.65	0.72	0.47	1.03
Single Ventricle w/ Arch Obstruction	1.72	0.77	3.71	6.20	1.94	21.52	0.84	0.02	13.07
Two Ventricles w/ Arch Obstruction	0.85	0.44	1.64	3.16	1.06	9.91	0.63	0.09	3.92
Age	0.98	0.92	1.05	0.97	0.87	1.07	1.21	1.02	1.45
Female	0.64	0.36	1.11	0.54	0.20	1.40	2.47	0.46	14.76